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1 Introduction 

1.1 Problem Description 

Recent developments ia audio and video systems require the use of multiple channel pro- 
cessing and reproduction with acoustic echo cancellers {AEC) and noise cancellers. In 
mioi-group "ideo conferencing systems multiple channel transmission leads to a better 
"localization* of the diverse people in the rooms. This enhances the intelligibility and 
naturalness of the speech. Stereo echo cancellation is needed in voice controlled stereo au- 
dio en video equipment. A multiple channel AEC can in general not be created by simple 
combination of multiple single AEC's, For two channel (stereo) A EC's this problem is 
well known and gets much attention In literature [11]. 

When the mmber of input signals In a multiple input acoustic echo canceller ia larger then 
the number of independent signal (and noise) sources in a multiple input acoustic echo 
canceller there is no longer a unique solution for the adaptive filters. This is the so called 
"non-uniqueness" problem. In practice however, the number of independent sources is 
always larger than the number of microphones due to noise and other disturbances 1 . When 
however the power- is not equally distributed over the independent sources, the problem 
might get badly conditioned, which is the main reason for the decrease in performance in 
the multiple input adaptive filters. In literature methods have been proposed to increase 
the performance by using Recursive Least Squares (RLS)-like algorithms (that have a 
huge computational complexity, even the most efficient implement at ions) and by adding 
non-linearities to the input signals (which might lead to audible artefacts in the output 
signal(s)). 

Similar problems arise in the multiple channel noise cancellation where noise that is picked 
up by a mk rophone is reduced using extra microphone signals as noise reference signal 
[7, 8]. The reference signals are filtered and subtracted from the (delayed) primary micro- 
phone signa.. 

i "Non-unic ueness" is the name used in literature to describe the problem fot the stereo acoustic echo 
canceller. Wb!n the input signals of a multiple input filter are filtered results of the same (single) source, 
there is no uiique solution for the adaptive filter [11], In practice this asms "non-uniqueness" ia not 
correct. As tte 'alternative' solution ia the stereo echo canceller case contains the product of a transfer 
function of a room and the inverse of another transfer function in the same room, it will in practice never 
be finite, mealing that there is in fact a unique solution. The presence of only a very small amount of 
noise (on for example the microphones) or the presence of at least as many sources as input signals to the 
multiple chazDiel adaptive filter are other re&aons why the "non-uniqueness" problem is in fact a problem 
of bad conditioning, see also [15]. 
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Here w* will present a multiple input frequency domain adaptive filter that make* use of 
the transformed inverse auto- and cross-correlation matrices to improve the performance 
of the adaptive filter* without a huge increase in computational complexity. 

1.2 Notation 

Signals are denoted by lower case characters, constants by upper case. Underlining is used 
for vectors, low case for time domain, and upper case for frequency domain. Matrices are 
denoted by caligraphic upper case, like X, or bold face upper case, like I. The dimension 
is put in superscript (e.g. ihsBxQ matrix X is given by for a square matrix the 
second dimension is omitted). Diagonal matrix are denoted by a double underline, like % 
with its diagonal denoted as E =» diag{£}. A subscript t, like m*, denotes the »'th version. 
The jsta element of ia is given by (is)«- Finally, appending [k] denotes the time index, 
(•)* denotes the transpose, (•)* the complex conjugate and (•)* the Hermitain transpose 
(complete conjugate transpose). 

2 Multiple Input Adaptive FIR Filter 
2.1 Adaptive Filter Structure 

A genertd multiple input adaptive FIR filter, depicted in figure 1, uses the S signals ar 0 [*] 
until x s _i[Jfe] to remove unwanted components related to these signals in the signal e p [*]. 
The signals xo[k] until *s_i{*J are input to S FIR filters Wo until Ws-x, with outputs 
e 0 [*3 un;il e s -i[fc]. The goal of the update algorithm is to adapt the coefficients of the 
FIR filters in such a way that the correlation between e p [fcj and the input signals *o(*i 
until zs.-itfcl is removed. For S > a > 0 the FIR filter W a performs the convolution of 
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Figure 1: Multiple input adaptive FIR filter. 
the signal and the coefficients w a ,D. . . w a ,/v-i of that filter. The output signal e»(fc] 
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of such a filter can be described as follows 



y-i 

4=s<J 



(i) 



with for S > a > 0 



*.[*-!] 



(2) 



(3) 



The output: of the multiple input adaptive filter is given by 



(4) 



These filter parts of the separate (adaptive) filters Wo until Vfs-x can be implemented 
efficiently in frequency domain with help of partitioning, block processing and FFTs as 
described in [1, 2, 3, 6, 10], 

Fox the ujxlate of the filter coefficients one can use S separate update algorithms, lb 
improve convergence behaviour the input signals can be decorrelated separately in the 
time donuiin by using RLS like algorithms, leading to a huge computational complexity. 
Complexii y reduction can be obtained by implementation in the frequency domain, with 
(Decoupled or Non-Uniform) Partitioned Block Frequency Domain Adaptive Filters (D 
or NU)*P:3PDAF ([5, 6, I0 f 12]). When there is correlation between the input signals of 
the filters this might however still lead to very bad convergence behaviour, see also [15], 
appendix C. 

2.2 Multiple Input FDAF 

We propose to use a partitioned algorithm in the frequency domain that reduces the effect 
of the crouscorrelafcion between the input signals on the algorithm's convergence behaviour. 
To reduce complexity we will use block processing with block length A to compute the sum 
of A consecutive updates with each iteration. We partition for S > a > 0 the coefficient 
vectors %?[IA] into g u parts of length Z with 3 g v =*= [N/Z] such that for S > j > 0 



(5) 



'When NjZ is not an integer, all y£ are extended to a length \N/Z\ - Z = g u *Z vector by appending 
\SJZ\ • Z - JY coefficients that are kept aero. 
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with fo) • S > 3 > 0 and Su > * > 0 



i 



(6) 



As we -vill use a. Fourier transform length L, with L > Z + A - 1, we define the input 
signal Fourier transforms for S > a > 0 as 

( x a [tA-L + l] \ 

* 

a? a [lA - 1] 

The diagonal matrices x£ [iA] contain, the vector X%[t A] as main diagonal, so for S > a > 0 

Jf PA] « diag{^[/A]>. (8) 

As we use an overlap-save method to compute the correlation involved in the adaptation 
proces in the frequency domain, the frequency domain transform of the residual signal 
vector equals 



tfpA]-** 



/ ft*-* ^ 
V L A U) J 



(9) 



The set of update equation* for the filter coefficients in the MFDAF algorithm can now 
be defined for g u > * > 0 by (for a full derivation see [15]) 



+ 1)A} 



+ 



with 



/ CSS-iPA])* \ 



(10) 



= 2a{-P s - L [lA])- 1 



and the tranaformation matrix Q s :^ s - L is given by 

/ g2,L qZL . . . qZ,l \ 

gS Z,S L _ 



(11) 



t 



qZ,1 gZ,L j 



(12) 
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= ( J* 0^-* ) • (^J -1 . (13) 

In order to be able to use a delay line for the diverse input signal transforms, we have to 
choose Z ia such a way that ZfA is integer, because then for all $ > i > 0 (£ L [lA-i2])* = 
QC*fC " *f ) j4 B* (see subsection 4.2.4 of [10]). 
The input channel's power matrix P SL [IA] is defined by 



P st [lA] m j£ {(X^llA})' (X^llAiy} 



(14) 



where 



«{Qtf[M])-a5*. l [Lt]} £i(X^[lA])'^[U)} 



(15) 



The expectation operator £{} of equation (14) has to be replaced by an estimation routine, 
which is th? subject of section 2.3. As we need the inverse oiP S L [lA] a lot of divisions are 
required (in an algorithm to compute the inverse of this sparse matrix). We can reduce the 
number of operations by using a direct inverse matrix estimation algorithm, as is discussed 
in subsection 2.4. Note that the inverse of T? a L [IA\ is also a sparse matrix with the same 
structure aid we define 



2a(P 5t [M])" 1 = 



(16) 



where for D < i < $ and 0 < j < S 

^m^d^iT^iA}}. ( l7) 

In figures 2 and 3 the algorithm is depicted. 

2.3 Power Estimation 
2.3.1 Est imatlon of 2a(P)~ 1 

Estimation of 2a(P) -1 can be done by estimating and inverting the resulting matrix. 
The special structure of this matrix (and its inverse) implies that we need § - 1 inversions 
of complex S x S matrices and two inversions of real S xS matrices. Because of estima- 
tion errors we need to limit the eigenvalues of and we obtain from [15] the following 
estimation routine 

HLU) - (i~ 7 )-pf;iUa-i)A] ( i8) 
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Figure 2: Calculate £ in MPBFDAF. 
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Pigure 3: Calculate w?[lA] in MPBFDAF. 
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where P,« in is the eigenvalueshift (aJid thus the 8 minimum for the eigenvalues of FfjU'^l- 
In appendix A the initialization for the estimation routine ia described. ' 

2.4 estimation of Inverse 
2.4.1 Plain Borders 

lb avcid the matrix inversion* we can also use a direct estimate of the inverse power 
matrix. To guarantee stability we have to use an algorithm that keeps the eigenvalues 
positive. In [15J, such an algorithm is derived and the results are summarised in table 1. 
In appendix A the initialization of the inverse power matrix is discussed. 



Initialisation 

1. S v ^2^» 3 ^ 

2. U 1 = V Praia ' 

for s = 1 to log 2 (5 u )do 

begin 

u - { u*- 1 -u*" 1 ) 

end 

3. Initialise power matrix (see appendix A). 
Iteration 



f US*-3,l mod 5. m l L ^ 
V Usu-U mod S u • * L J 



4. Calculate Qft 6 ^])" 1 

-Q***[M] • (D L [M))-* ■ (Q^M)* 

Tfcble 1: Direct inverse power update with limits. 



2.4.2 Power Matrix Dimension Reduction 

Partitianing of the coefficient vectors implies a reduction in the dimension of the power 
vector and an increase in both the adaptive filter update rate and the update rate of 
the power matrix estimate (assuming that the block length is in the same order as the 



8 



9 

partition :ength) s . J Due to the increased rate and the reduced dimension we can make a 
much moie accurate estimate of the inverse power matrix. Besides that, the tracking of 
the non-si ationarities in the input signal power is improved and, due to the increased rate 
in the up iate part, alao the tracking of the non-etationarities in the unknown path, see 
also 2.5.4 and 4.4,3 in [10]. 

2*4.3 E ^coupling Power Estimate 

Partitioning of the coefficient vector Increases the computational complexity. When the 
main cauue of trouble is the tracking of the input signal power, we can reduce the dimen- 
sion of tbe input signal power matrix and increase ita update rate independent from the 
partition length of the coefficient vector. We then can obtain part of the advantages of 
partitioning with a much smaller increase in complexity, see subsection 6,4.2 in [10]* 

3 Computational Complexity 

To obtain a rough estimate of the effect on the computational complexity of the proposed 
algorithm, we will compare our MFDAF algorithm with using S separate DPBFDAF algo- 
rithms. We define (?« and G+ as tbe number of operations needed for a multiplication and 
an addition. Gp and <?(p)-i 1 denote the number of operations needed for the estimation 
of P, indirect (subsection 2.3.1) or direct (subsection 2.4). F{L} denotes the number of 
operations needed for a length L real input PPT. The total number of real operations for 
the update part per block of A samples is summarised in table 2 (see [15]). 





Separate DPBFDAFs 


Multiple Input BFDAF 


;i)PFTs 


((! + *„) -S+D-W} 


p-l- ' 

indirect 
[1!], 88. 3.1.3 




Gr{S}+G-i{S} 


direct 
[1!.], as. 3.2.4 










£-5'((2.5-l) + 1) -G+) 













Table 2; Aeai operations per A samples for update algorithms. 

When \/e compare the Multiple Input BFDAF to the separated channels case, we see 
that thisre is a difference in the number of elementwise multiplications and in the power 
update. The main computational load however, is caused by the computation of the FFTs, 
so this ncrease will only make a small relative difference (depending on S) For example, 
when *e use split-radix FFTs whose complexity is given in table 3 y that » G* = 1 
(see [1C]) and that G, = 5, and take N = 2048, A = Z = 512, L = 1024, g u = 4 and 

3 The rate depeada only on the block length, but this block length is strongly related to the partition 
length ii order to lceep the computational complexity aft low as possible. We can of course also reduce 
the blocl; length. and thus increase the rate without partitioning. This however leads to a huge increase in 
computational complexity (see (10]). 
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r 
jj 


R. 1 in G*) 


" Fx ( in GU) 


2 


o 


2 


4 


o 


6 


g 


2 


20 


16 


12 


58 


32 


42 


154 


64 


122 


392 


128 


328 


950 


256 


822 


2240 


512 


1984 


5158 


1024 


4646 


11680 


2048 


10656 


26086 


4096 


24038 


57632 



Table 3: Operations ibr aplit-r&dix real input FFT. 



the caws S = 2 and 5 - 3, we get with help of table 2, the results of table 4. We 
can se; that there is an increase of approximately 10 % per extra input channel of the 
multiple input BPDAF compared to the separate DPBFDAF's. As the update part takes 
appro* imately 50 % of the computational toad in the BFDAF case (the other half is needed 
for the actual filtering operation), the increase in computational complexity for the whole 
multiple channel adaptive filter is about 5 % per extra input channel. For 5 > 4 the 
computational complexity of the direct inverse power estimate is much smaller then the 
computational complexity of the power estimate and matrix inversion method. 



Operation 


Sej 


sarate DPBFDAFs 


Mul 


Ltiple Input BFD 


AF 




S*>2 


5 = 3 


5 = 4 


5 = 8 


5 = 2 


5 = 3 


5 = 4 


5 = 8 


(I)FFTs 


351 


510 


670 


1307 


351 


510 


670 


1307 


P" 1 (indirect) 


30 


45 


60 


120 


50 


146 


904 


6712 


T- 1 (direct) 


52 


78 


104 


208 


81 


170 


293 


1125 


— P^ X 


4 


6 


8 


"' 16 


16 


42 


80 


352 




48 


72 


96 


192 


48 


72 


96 


" 192 




8 1 


12 


16 


32 


8 


12 


16 


32 


TTctal (indirect) 


441 


645 


850 


1667 


473 


782 


1766 


8595 


Total (direct) 


463 


678 


894 


1755 


504 


806 


1155 


3008 


Optimum 


441 


645 


850 


1667 


473 


782 


1155 


3008 


Increase 
for MFDFAF 










7% 


21% 


36% 


80% 



Table 4: Example real operations per sample for update algorithms. 
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4 Applications 

4*1 St ereo Acoustic Echo Canceller 

A full stereo communication requires four stereo AECs, two on the near end side and two 
on the fix end side. In figure 4 one of those echo cancellers is depicted. Note that on 
each sid« we can combine the input signal delay-lines, the FFTs and the multiplication 
by the inverse power matrix of the two echo cancellers, which implies that the relative 
extra computational complexity for removing the crosscorretation is even furter reduced. 
The performance of the AECs is further improved by adding Dynamic Echo Suppressors 



"far'fdofn 



ID- 



ID- 



Adaptive 
Filter 




Figure 4: Isolated Adaptive Filter in Stereo Echo Cancellation. 
(DBS) aid programmable filters, see [4, 13, 14], leading to the configuration in figure 5. 
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Figure S: Isolated Adaptive Filter with Dynamic Echo Suppressor. 
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4.2 Multiple Input Adaptive Noise Canceller 



A second application for the multiple input adaptive filter is the multiple input noise 
canceller. The first configuration ia depicted in figure 6. The microphones are placed in 
a cross 4 . We use a selection mechanism to select the microphone were the desired signal 
(spoaker) ia located and use the other microphone signals (xo[*l until x *\ k \ ( or i mxe 
general, x S -i[k\)) to cancel the other signal sources. 




Figure 6: Adaptive ATofae Canceller with a Cross ConSgw&tion. 

A second configuration tries to select the desired signal from an array of microphones using 
a dehy-sum beamformer, and creates noise (unwanted signals) references by subtraction, 
see figure 7. Also in the multiple input Noise Canceller we can apply Dynamic Echo 
Suppressors (which is in feet not surpressing an echo, but is similar to the DBS in the 
AECs case) and programmable filters to improve performance (which is not depicted 
here). An extra problem is that the inputs of the filters may contain some elements of the 
desired signal ("signal leakage", because the delay-sum beamformer is not perfect and/or 
the directional microphones in the cross configuration are not ideal). When the desired 
signal is a speech signal, a speech detector can be used to improve the behaviour of the 
MFDj^F. Another solution to suppress the influence of the signal leakage ia given in [9]. 



4 Th: number of microphones is of course not restricted to font. 
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Figure 7: Adaptive Noise Canceller with a Linear Array. 

5 Conclusions 

A Multiple input Frequency Domain Adaptive Filter (MFDAF) algorithm is presented 
that uses partitioning and decorrelation techniques in the frequency domain to improve 
convergence behaviour, especially in highly correlated and/or non-stationary environment. 
The condit ioning of the problem is greatly improved by using estimates of the transformed 
auto- and crosscorreiation matrices (see [16], appendix C), Efficient algorithms for the 
estimation of the inverse transformed auto- and crosscorreiation matrices are presented so 
using MFDAF does not lead to a huge increase in computational complexity compared to 
using multiple single input adaptive filters. The main applications are the stereo Acoustic 
Echo Canceller and the multiple input adaptive Noise Canceller. 

A Initialization of power matrix 
A.l Estimation of P 

Equation (18) that describes the power matrix update, was given by 

H%»*\ = (1-7)'FS4.[(«-1M]+ (19) 

y ' (djl ' (x5 " (u)r ' < x * ^ W + if ' • 

Prom th< above equation we see that (for i £ /) {X $LA [{1 - i)A])* • (X s ^{(l - i)yl])* 
is weighted by a factor (1 - 7 )< compared to (X SmL *[lA]) m • (X* L > L [tA])K As we do not 
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have any a priori knowledge of the input signals, we cannot make an accurate guess for 
the ijtifcial power matrix Pftfml-A]. When WW assume that it equal the all *ero matrix, 
Hie matrix Pj&IM] is a factor i - (1 - too small (in a stationary environment). 
We therefore propose to use Pf^mjaitMli with 

pa^iMi - ^^L^i • p &[M). (20) 

This implies that 

As this takes a few extra multiplications and divisions, we propose to do this until 
1 — (1 — 7)'+ l is close to 1 (e.g. > 0.S5), and then switch to the not corrected formula 
(with the corrected power matrix as first guess). 

A.2 Intimation of Inverse 

Fos the direct inverse power matrix estimation we use in table 1 



with 



and 



(22) 
(23) 
(24) 



Here we d 3 not have the possibility of an initial guess with the «ero matrix for (Pf A))" , 
because tlien the algorithm would always yield a zero matrix as result. We therefore sug- 
gest to up; a momentaneous estimate of the diagonal of Pf [am[0], so 



(P»jttm,iiil["~4l)""* 



2al 



\ 



0* 



0* 



(25) 
-l 



2a 



which meais that we have to calculate S ■ L divisions. To get a fair weighting of the 
input vectors in the first iterations, also here we have to adjust the multiplication Cactors 
in the update algorithm. In the first iterations 7 must be chosen larger to get a better 
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estimate of the full power matrix, so we propose to use TinU ■ - 7) 1 ™tead of 7 as long 
as Tin* ■ ( 1 — 7)' » larger than 7. 3 



(Hfa.mtW = 1 - Tu|it 1 (1 - 7) * ■ CPgftgtK* " ^P"* " 

QW[1A] • (fi^Mj)- 1 • <Q*^[M])* (26) 



with 



l-Hnit-a-T)' 



(27) 



and 



(28) 



Tintt • (1 - 7)' 



B Instabilities 

The power matrix estimation and inversion algorithms do not give an absolute guarantee 
that the bop of the adaptive filter is stable. The estimation errors and approximations 
used to mike the correlation matrices circulant may result in a temporarily feed bad? g*"P 
for one of the eigenxnodii in the adaptive filter that is to large, especially when the input 
signal is highly non-stationary. To guarantee stable operation of the MFDAP algorithm 
(including programmable filters) we therefore need to check if such instabilities occur and 
reinitialize the power matrix (and filter coefficients) if necessary* The coefficients calculated 
by the adaptive filter are only tranafered to the programmable filters (that calculate the 
actual output of the system), when this leads to a reduction in the power of the output 
signal (residual signal) (see [13]). Instability therfore means that the output generated 
by the adaptive filter explodes, while the actual output (generated by the programmable 
filters) dots not do that. 

We detect Instabilities by comparing the mean squared value of the adaptive filter output 
by the mam squared value of the programmable filters output, When this ratio exceeds a 
certain thieshold (»1) and there is relevant input 6 the adaptive filter is reset (the power 
matrix is reinitialised and the coefficients are put to aero). 



Note that y mi% cannot be 1, as we divide by 1 - ^ m the first iteration! 

This is cliecked by demanding that the ratio of the power of the adaptive filter output and the power 
of the signal ;hat must be corrected exceeds the same threshold. When this happens, the signal to correct 
is much smaller than the signal it is corrected by, this certainly does not improve the result. 
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1. 



Stereo acoustic echo canceller comprising decorrelatian means. 



2. Adaptive noise canceller comprising decorrelation means. 

5 3. Decorrelation means for use in a stereo acoustic echo canceller as claimed in 

claim 1 



4. 



Decorrelatioa means for use in a adaptive noise canceller as claimed in claim 2. 



10 5. 



Method of decorrelation of at least two signals. 
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